
#ifndef ROPUFU_MULTI_SENSOR_MONTE_CARLO_DUMPER_HPP_INCLUDED
#define ROPUFU_MULTI_SENSOR_MONTE_CARLO_DUMPER_HPP_INCLUDED

#include "../distributed/matrix.hpp"
#include "../format/matstream_v4.hpp"

#include <aftermath/algebra/mixed_numeral.hpp>
#include <aftermath/probability/empirical_measure.hpp>

#include <cstdint>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>

namespace ropufu
{
    namespace monte_carlo
    {
        namespace detail
        {
        }

        template <typename t_stopping_time_type, typename t_mixed_numeral_type = aftermath::algebra::mixed_numeral<100>>
        struct dumper
        {
            typedef t_stopping_time_type stopping_time_type;
            typedef t_mixed_numeral_type mixed_numeral_type;
            typedef dumper<stopping_time_type, mixed_numeral_type> type;

            //const std::size_t count_streams = statistic_type::count_streams;

            typedef aftermath::probability::empirical_measure<stopping_time_type, double> empirical_time_type;
            typedef aftermath::probability::empirical_measure<double, double> empirical_double_type;
            typedef aftermath::probability::empirical_measure<bool> empirical_bool_type;

        private:
            double m_mu, m_theta, m_eta;

        public:
            explicit dumper(double mu, double theta, double eta)
                : m_mu(mu), m_theta(theta), m_eta(eta)
            {
            }

            static bool does_exist(const std::string& filename)
            {
                bool result = true;
                std::ifstream filestream;
                filestream.open(filename.c_str(), std::ios::in | std::ios::binary);
                if (filestream.fail()) result = false;
                else filestream.close();
                return result;
            }

            template <typename t_statistic_type>
            std::string build_filename(const t_statistic_type& procedure,
                std::uint_fast64_t n, std::size_t n_affected, const std::vector<double>& threshold,
                bool do_stop_on_reset) const
            {
                double expected_shift = this->m_theta - this->m_mu;
                double actual_shift = this->m_eta - this->m_mu;

                std::ostringstream filename;

                // Warning: if the path, along with the ".../theta_1_00/", or ".../theta_0_50/" subfolder does not exist, the program (dumping) might not work properly.
#ifdef _WIN32
                filename << "E:/Numerical/MATLAB Data/"; // TODO: move path selection elsewhere.
#else
                filename << "/local_hd/gsokolov/Numerical/";
#endif
                filename << "xtheta_" << mixed_numeral_type(expected_shift) << "/";

                filename << "ms";
                filename << "_k=" << t_statistic_type::count_streams;
                filename << "_n=" << (n / 1000) << "e3_" << procedure;
                if (do_stop_on_reset) filename << "_repeated";
                filename << "_a_from_" << mixed_numeral_type(threshold.front()) << "_to_" << mixed_numeral_type(threshold.back());
                filename << "_" << "lin";
                filename << "_in_" << threshold.size();

                if (actual_shift == 0.0 || n_affected == 0) filename << "_arl";
                else if (actual_shift == expected_shift) filename << "_" << n_affected << "add";
                else filename << "_" << n_affected << "xxx";

                filename << ".mat";

                return filename.str();
            }

            void dump(const std::string& filename, std::uint_fast64_t n, std::size_t n_affected,
                const std::vector<double>& threshold, const std::vector<empirical_time_type>& run_length) const
            {
                std::vector<empirical_double_type> aux;
                this->dump(filename, n, n_affected, threshold, run_length, aux, false);
            }

            void dump(const std::string& filename, std::uint_fast64_t n, std::size_t n_affected,
                const std::vector<double>& threshold, const std::vector<empirical_time_type>& run_length,
                const std::vector<empirical_double_type>& aux, bool do_aux = true) const
            {
                //double expected_shift = this->m_theta - this->m_mu;
                double actual_shift = this->m_eta - this->m_mu;

                // Dump results into a MATLAB MAT-v4 file.
                distributed::matstream<4> out(filename);
                out.clear();

                // Name of the vector. Should contain only numbers and letters, and start with a letter.
                out << "A";
                // Convert the vector of thresholds to a matrix, and write it to the file.
                out << distributed::matrix_helper::create_row_vector(threshold);


                if (actual_shift == 0.0 || n_affected == 0)
                {
                    out << "T";
                    // Unfold the empirical measure; not (hard-drive-wise) memory efficient, but straightforward to read and store.
                    distributed::matrix<stopping_time_type> observed_run_length(n, threshold.size());
                    for (std::size_t j = 0; j < threshold.size(); j++)
                    {
                        std::size_t flat_pos = 0;
                        if (run_length[j].count_observations() != n) std::cerr << "Size error.";
                        // Iterate over all observed stopping times for current threshold; <item> is a pair (observed value, # of times it has been observed).
                        for (const auto& item : run_length[j].get_data())
                        {
                            // If current item is, e.g., (tx, 15), write (tx, tx, \cdots , tx)', a 15-by-1 vector to the right position in the matrix.
                            for (std::size_t i = 0; i < item.second; i++) observed_run_length.write(flat_pos++, j, item.first);
                        }
                    }
                    // Write the matrix to the file.
                    out << observed_run_length;
                }
                else
                {
                    stopping_time_type min_stop = run_length[0].get_min();
                    stopping_time_type max_stop = run_length[0].get_max();
                    for (std::size_t j = 1; j < threshold.size(); j++)
                    {
                        if (run_length[j].get_min() < min_stop) min_stop = run_length[j].get_min();
                        if (run_length[j].get_max() > max_stop) max_stop = run_length[j].get_max();
                    }
                    std::size_t n_values = (max_stop + 1) - min_stop;
                    std::vector<stopping_time_type> observed_stopping_values(n_values);
                    for (stopping_time_type i = 0; i < static_cast<stopping_time_type>(n_values); i++) observed_stopping_values[i] = min_stop + i;

                    out << "T_values";
                    // Convert the vector of thresholds to a matrix, and write it to the file.
                    out << distributed::matrix_helper::create_column_vector(observed_stopping_values);

                    out << "T_counts";
                    distributed::matrix<stopping_time_type> observed_stopping_counts(n_values, threshold.size());
                    observed_stopping_counts.fill(0);
                    for (std::size_t j = 0; j < threshold.size(); j++)
                    {
                        const auto& run_length_for_threshold = run_length[j];
                        if (run_length_for_threshold.count_observations() != n) std::cerr << "Size error.";

                        // Iterate over all observed stopping times for current threshold; <item> is a pair (observed value, # of times it has been observed).
                        for (const auto& item : run_length_for_threshold.get_data())
                        {
                            std::size_t row_index = item.first - min_stop;
                            observed_stopping_counts.write(row_index, j, static_cast<stopping_time_type>(item.second));
                        }
                    }
                    // Write the matrix to the file.
                    out << observed_stopping_counts;
                }

                if (do_aux)
                {
                    out << "Aux";
                    // Unfold the empirical measure; not (hard-drive-wise) memory efficient, but straightforward to read and store.
                    distributed::matrix<double> observed_aux(n, threshold.size());
                    for (std::size_t j = 0; j < threshold.size(); j++)
                    {
                        std::size_t flat_pos = 0;
                        if (aux[j].count_observations() != n) std::cerr << "Size error.";
                        // Iterate over all observed stopping times for current threshold; <item> is a pair (observed value, # of times it has been observed).
                        for (const auto& item : aux[j].get_data())
                        {
                            // If current item is, e.g., (tx, 15), write (tx, tx, \cdots , tx)', a 15-by-1 vector to the right position in the matrix.
                            for (std::size_t i = 0; i < item.second; i++) observed_aux.write(flat_pos++, j, item.first);
                        }
                    }
                    // Write the matrix to the file.
                    out << observed_aux;
                }
            }
        };
    }
}

#endif // ROPUFU_MULTI_SENSOR_MONTE_CARLO_DUMPER_HPP_INCLUDED
